function cgm2

%  check convergence of CGM for solving Ax=b 
%  and compare different error measures in the process

% clear all previous variables and plots
clear *
clf

% pick exact solution
% ic=1: b=ones
% ic=2: b=random
% ic=3: b=fourier series sol
ic=3;

% loop for increasing N
for ie=1:3

	% set parameters and define matrix A
	N=25*2^(ie-1);
	M=N;
	nm=N*M
	lambda=1;
	A=matrix(lambda,N,M);

	% define b
	if ic==1
		sol=ones(nm,1);
	elseif ic==2
		sol=rand(nm,1);
	elseif ic==3
		sol=fourier(N,M)';
	end;
	b=A*sol;

	start_error=norm(sol,2);

	fprintf('\n n = %i    Initial Error = %e\n\n',nm, start_error)

	% use CGM to solve matrix equation '
	v=cgm(A,b,sol,ie);

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v=fourier(NN,M)
modes=100;
a=1; b=1;
k=b/(M+1);  h=a/(NN+1);
for n=1:modes
	a1=cos(3*n*pi/(4*a))-cos(n*pi/(4*a));
	an(n)=-2*a1/(n*pi*sinh(n*pi*b/a));
end;

for j=1:M
	y=j*k;
	for i=1:NN
		x=h*i;
		s=0;
		for n=1:modes
			s=s+an(n)*sinh(n*pi*y/a)*sin(n*pi*x/a);
		end;
		ll=(j-1)*NN+i;
		v(ll)=s;
	end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction for creating the matrix A
function A=matrix(lambda,N,M)

nm=N*M;

% generate the diagonal elements
D = sparse(1:nm,1:nm,2*(1+lambda^2)*ones(1,nm),nm,nm);

% generate the subdiagonal elements
SD = sparse(2:nm,1:nm-1,-ones(1,nm-1),nm,nm);

for i=N+1:N:nm-1
	SD(i,i-1)=0;
end;

% generate the sub-subdiagonal elements
SSD=sparse(N+1:nm,1:nm-N,-lambda^2*ones(1,nm-N),nm,nm);

% use symmetry of A to complete construction
A=D+SD+SD'+SSD+SSD';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction for the CGM
function x=cgm(A,b,sol,ie)

% set CGM parameters
tol=0.0000001;

nm=length(b);
x=zeros(nm,1);
tic
% start iteration
r=b-A*x;
d=r;
rr=dot(r,r);
error=1;
counter=0; 
err=1;
while err>tol
	counter=counter+1; iter(counter)=counter;
	if counter==1
		beta=0;
	else
		beta=rr/rr0;
	end;
	d=r+beta*d;
	q=A*d;
	alpha=rr/dot(d,q);
	x=x+alpha*d;
	r0=r;
	r=r-alpha*q;
	rr0=rr;
	rr=dot(r,r);
	error2(counter)=norm(alpha*d);
	error3(counter)=norm(r);
	error(counter)=norm(x-sol); err=error(counter);

	if counter>2
		if (error2(counter)<1e-6) & (error2(counter-1)>1e-6)
			fprintf('\n  Number of CGM Iterations = %i     Error =  %e    Time = %e \n\n',counter,error2(counter),toc)
		end;
	end;

end;

% plot results

% get(gcf)
%set(gcf,'Position', [377 468 558 548]);
%set(gcf,'Position', [1100 368 589 505]);
plotsize(558, 548)

subplot(3,1,ie)

loglog(iter,error,'k')
hold on
loglog(iter,error2,'r')
loglog(iter,error3,'b')

% axes and legend used in plot
ylabel('Error','FontSize',14,'FontWeight','bold')
if ie==3
	xlabel('Iteration Steps','FontSize',14,'FontWeight','bold')
	text(2,0.00001,'n = 10000','FontSize',14,'FontWeight','bold')
elseif ie==1
	legend(' Error',' Iteration Error',' Residual',1);	
	text(2,0.00001,'n = 625','FontSize',14,'FontWeight','bold')
		say=['Solving Laplaces Equation Using the CGM'];
	title(say,'FontSize',14,'FontWeight','bold')
elseif ie==2
	text(2,0.00001,'n = 2500','FontSize',14,'FontWeight','bold')
end;

% have MATLAB use certain plot options (all are optional)
grid on
set(gca,'MinorGridLineStyle','none')
axis([ 1 1000 0.00000001 100])
set(gca,'YTick',[0.00000001 0.000001 0.0001 0.01 1.0 100])
set(gca,'FontSize',14);
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 

hold off

%fprintf('\n  Number of CGM Iterations = %i     Error =  %e \n\n',counter,error(counter))
%cond(full(A),2)

% subfunction plotsize    '
function plotsize(width,height)
siz=get(0,'ScreenSize');
bottom=max(siz(4)-height-95,1);
set(gcf,'Position', [2 bottom width height]);